基于 R 语言的科研论文绘图技巧详解（5）

点击下方公众号，回复资料分享，收获惊喜

简介

在查阅文献的过程中，看到了几幅非常不错的出版图，今天就跟着小编一起学习下，他们是怎么使用 R 绘制出来的。

今天主要介绍 第五幅图（E) —— 带置信区间的拟合曲线图及带误差项的散点图。这个图在科研绘图中较为常用，例如：绘制产品失效时间的散点图。

前四幅图的详细代码介绍可见：基于 R 语言的科研论文绘图技巧详解（4）、基于 R 语言的科研论文绘图技巧详解（3）基于 R 语言的科研论文绘图技巧详解（2）基于 R 语言的科研论文绘图技巧详解（1） 。最后一幅图会随后继续介绍，读者在学习过程中，可以将内部学到的知识点应用到自己的图形绘制中。推文已经将主要知识点进行罗列，更有利于读者学习和查阅。

那我们来看看，作者是怎么实现这个功能的吧，本文知识点较多，大家耐心学习，建议自己实践。对应代码、数据可在 GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R^[1] 中找到。

主要知识点

• 绘制散点图、丝带形状图；
• 绘制横向、纵向误差图；
• 学会小技巧：展示轴外部的图形。

绘图

加载包

首先加载一些需要使用到的包。

library(ggplot2) # Grammar of graphics

设置主题

接下来，为了方便起见，作者在绘图前设置好了主题，并将该函数命名为 my_theme。

这一部分在第一篇推文给出，代码将在文末中完整代码给出。

手动修改大部分面板，具体可以参考本篇文章^[2]。或者观看我在 B 站发布的《R 语言可视化教程》，里面也有一些简单主题设置介绍。

导入数据

首先使用 read.csv() 导入三个基因相关数据，最后合并到 data_E 数据框中。前 6 行数据如下所示。

data_Ea = read.csv("./data_Ea.csv")
data_Eb = read.csv("./data_Eb.csv")
data_Ec = read.csv("./data_Ec.csv")
data_E = data.frame("gene"=c(rep("gene a",nrow(data_Ea)),
                             rep("gene b",nrow(data_Eb)),
                             rep("gene c",nrow(data_Ec))),
                    "t"=c(data_Ea$t,data_Eb$t,data_Ec$t),
                    "C"=c(data_Ea$C,data_Eb$C,data_Ec$C),
                    "pos_err"=c(data_Ea$err,data_Eb$pos_err,data_Ec$pos_err_C),
                    "neg_err"=c(data_Ea$err,data_Eb$neg_err,data_Ec$neg_err_C),
                    "pos_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$pos_err_t),
                    "neg_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$neg_err_t)
)

head(data_E)
#    gene  t    C pos_err neg_err pos_err_t neg_err_t
# 1 gene a 12 0.27    0.03    0.03        NA        NA
# 2 gene a 24 0.19    0.02    0.02        NA        NA
# 3 gene a 36 0.17    0.01    0.01        NA        NA
# 4 gene a 48 0.13    0.04    0.04        NA        NA
# 5 gene a 60 0.09    0.03    0.03        NA        NA
# 6 gene a 72 0.10    0.01    0.01        NA        NA

gene 是分类变量，t 表示时间将作为 x 轴。C 表示正常水平下的结果（小编对基因方面不是很懂，可能理解的不对）。pos_err 和 neg_err 是正向与反向误差。

自定义函数

首先，定义三个函数，以便稍后拟合散点图。这里的参数具体怎么得到的，作者没有介绍，小编不是很了解。或许是使用了某种参数估计的方法吧。

f1 = function(t) 0.2*exp(-t/48)
f2 = function(t) 0.3*exp(-t/60)
f3 = function(t) 0.4*exp(-t/72)
t = seq(0,96,1)

# 构造 f2 函数的数据
ribbon = data.frame(
  "f2" = 0.3*exp(-t/60),
  "t" = t
)

绘制图形

定义图形形状分别为实心正方形、圆形、三角形 manual_pch =c(15,16,17) 。之后简单绘制下散点图geom_point()和 丝带形状的图 geom_ribbon()。

注意：这里的 factor() 将 gene 变成离散形式。由于 geom_ribbon() 中使用的数据不同，所以要重新指定下新数据。

manual_pch =c(15,16,17) 定义图形形状
ggplot(data=data_E) +
  geom_point(aes(x=t,y=C,pch=factor(gene)),size=2) +
  geom_ribbon(data=ribbon, aes(x=t,ymin=0.85*f2,ymax=1.15*f2),
              fill="black",alpha=0.1) 

在这个基础上使用 stat_function() 添加三个函数，两条虚线，一条实线。

stat_function(fun=f1, geom="line",linetype="dashed") +
  stat_function(fun=f2, geom="line") +
  stat_function(fun=f3, geom="line",linetype="dotted")

使用 geom_errorbar() 和 geom_errorbarh() 添加误差棒（纵向与横向）。注意内部参数不同：（ymin，ymax） 和 （xmin，xmax）。

geom_errorbar(aes(x=t,ymin=C-neg_err, ymax=C+pos_err), width=2) +
geom_errorbarh(aes(y=C,xmin=t-neg_err_t,xmax=t+pos_err_t))

最后就是主题的设置了，包括运用了自定义的主题 my_theme()。

scale_shape_manual(values=manual_pch) +
  my_theme() + theme(legend.title=element_blank())+
  theme(legend.position = c(0.9,0.95)) +
  scale_x_continuous(expand = c(0, 0), 
                     breaks = c(seq(0,96,12)),
                     limits = c(0,96)
  ) +
  scale_y_continuous(expand = c(0, 0),
                     breaks = c(seq(0,0.4,0.05)),
                     limits = c(0,0.4)
  ) +
  theme(
    panel.grid.major = element_line("gray95", size = 0.1),
    # putting label closer to axis bc exponent makes it bigger 
    axis.title.y = element_text(margin = margin(r = -9)) 
  ) +
  xlab("time (h)") +
  ylab(expression(paste("concentration (mmol",~cm^-3,")"))) 
 

但是，最右边几个点并没有显示清楚，这里作者运用了一个小技巧，将这些点清晰的呈现出来。来看看这个代码的效果：

 coord_cartesian(clip = "off") # to allow for plotting outside axes

完整代码

# Panel E ----
library(ggplot2) 
base_size = 12 
my_theme <-  function() {
  theme(
    aspect.ratio = 1,
    axis.line =element_line(colour = "black"),  
    
    # shift axis text closer to axis bc ticks are facing inwards
    axis.text.x = element_text(size = base_size*0.8, color = "black", 
                               lineheight = 0.9,
                               margin=unit(c(0.3,0.3,0.3,0.3), "cm")), 
    axis.text.y = element_text(size = base_size*0.8, color = "black", 
                               lineheight = 0.9,
                               margin=unit(c(0.3,0.3,0.3,0.3), "cm")),  
    
    axis.ticks = element_line(color = "black", size  =  0.2),  
    axis.title.x = element_text(size = base_size, 
                                color = "black", 
                                margin = margin(t = -5)), 
    # t (top), r (right), b (bottom), l (left)
    axis.title.y = element_text(size = base_size, 
                                color = "black", angle = 90,
                                margin = margin(r = -5)),  
    axis.ticks.length = unit(-0.3, "lines"),  
    legend.background = element_rect(color = NA, 
                                     fill = NA), 
    legend.key = element_rect(color = "black",  
                              fill = "white"),  
    legend.key.size = unit(0.5, "lines"),  
    legend.key.height =NULL,  
    legend.key.width = NULL,      
    legend.text = element_text(size = 0.6*base_size, 
                               color = "black"),  
    legend.title = element_text(size = 0.6*base_size, 
                                face = "bold", 
                                hjust = 0, 
                                color = "black"),  
    legend.text.align = NULL,  
    legend.title.align = NULL,  
    legend.direction = "vertical",  
    legend.box = NULL, 
    panel.background = element_rect(fill = "white", 
                                    color  =  NA),  
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.title = element_text(size = base_size, 
                              color = "black"), 
  ) 
  
}

data_Ea = read.csv("./data_Ea.csv")
data_Eb = read.csv("./data_Eb.csv")
data_Ec = read.csv("./data_Ec.csv")
data_E = data.frame("gene"=c(rep("gene a",nrow(data_Ea)),
                             rep("gene b",nrow(data_Eb)),
                             rep("gene c",nrow(data_Ec))),
                    "t"=c(data_Ea$t,data_Eb$t,data_Ec$t),
                    "C"=c(data_Ea$C,data_Eb$C,data_Ec$C),
                    "pos_err"=c(data_Ea$err,data_Eb$pos_err,data_Ec$pos_err_C),
                    "neg_err"=c(data_Ea$err,data_Eb$neg_err,data_Ec$neg_err_C),
                    "pos_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$pos_err_t),
                    "neg_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$neg_err_t)
)

head(data_E)

f1 = function(t) 0.2*exp(-t/48)
f2 = function(t) 0.3*exp(-t/60)
f3 = function(t) 0.4*exp(-t/72)
t = seq(0,96,1)

ribbon = data.frame(
  "f2" = 0.3*exp(-t/60),
  "t" = t
)
head(ribbon)

manual_pch =c(15,16,17) # available pch: type ?pch 

panel_E <- ggplot(data=data_E) +
  geom_point(aes(x=t,y=C,pch=factor(gene)),size=2) +
  geom_ribbon(data=ribbon, aes(x=t,ymin=0.85*f2,ymax=1.15*f2),
              fill="black",alpha=0.1) +
  stat_function(fun=f1, geom="line",linetype="dashed") +
  stat_function(fun=f2, geom="line") +
  stat_function(fun=f3, geom="line",linetype="dotted") +
  geom_errorbar(aes(x=t,ymin=C-neg_err, ymax=C+pos_err), 
                width=2) +
  geom_errorbarh(aes(y=C,xmin=t-neg_err_t,xmax=t+pos_err_t))+
  scale_shape_manual(values=manual_pch) +
  my_theme() + theme(legend.title=element_blank())+
  theme(legend.position = c(0.9,0.95)) +
  scale_x_continuous(expand = c(0, 0), 
                     breaks = c(seq(0,96,12)),
                     limits = c(0,96)
  ) +
  scale_y_continuous(expand = c(0, 0),
                     breaks = c(seq(0,0.4,0.05)),
                     limits = c(0,0.4)
  ) +
  theme(
    panel.grid.major = element_line("gray95", size = 0.1),
    # putting label closer to axis bc exponent makes it bigger 
    axis.title.y = element_text(margin = margin(r = -9)) 
  ) +
  xlab("time (h)") +
  ylab(expression(paste("concentration (mmol",~cm^-3,")"))) +  
  coord_cartesian(clip = "off") # to allow for plotting outside axes

panel_E

小编有话说

本文主要学到的知识点如下：

1. 使用 geom_point() 绘制散点图， geom_ribbon()绘制丝带形状图；
2. 使用 stat_function() 添加函数曲线；
3. 使用 geom_errorbar() 和 geom_errorbarh 添加误差棒（纵向与横向）；
4. 使用 coord_cartesian(clip = "off") 允许展示外轴的图形。

看完这篇文章，相信你的技能包又多了点东西。记得实操噢！看了不代表就会了～ 如果觉得内容有用的话，小编写的有心的话。给小编来杯咖啡吧！

参考资料

[1]

GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R: https://github.com/marco-meer/scifig_plot_examples_R

[2]

文章: https://ggplot2.tidyverse.org/reference/theme.html

推荐： 可以保存以下照片，在b站扫该二维码，或者b站搜索【庄闪闪】观看Rmarkdown系列的视频教程。Rmarkdown视频新增两节视频（写轮眼幻灯片制作）需要视频内的文档，可在公众号回复【rmarkdown】